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In a recent experiment [D. J. McCarron et al, Phys. Rev. A, 84, 011603(R) (2011)] a two- 
species 87 Rb- 133 Cs Bose-Einstein condensate was formed and three distinct regimes of density 
distributions observed depending on relative atom numbers of the two species. To investigate these 
theoretically, we obtain time-independent solutions of the trapped two-species condensate through 
zero-temperature mean-field simulations. We find the results to be sensitive to experimentally 
relevant shifts in the potentials in both longitudinal and transverse directions, and observe a range 
of structures, including 'ball and shell' formations and axially /radially separated states. We find 
good overall agreement with the experimental results for all regimes. 
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I. INTRODUCTION 

Since the successful realization of an atomic Bose- 
Einstein condensate (BEC) composed of two different 
hyperfine spin states of 87 Rb [1], experimental and the- 
oretical work has advanced greatly in the field of two- 
component BECs. These have been produced using dif- 
ferent atomic species [2-5] , different isotopes of the same 
atom [6], and a single isotope in two different hyper- 
fine spin states [1, 7-14]. Spinor condensates, which 
have at least three components with internal spin de- 
grees of freedom, are also generating much current in- 
terest (see [15] for a review). Due to the interactions 
between the species, a key feature of two-species BECs 
is their potential to exhibit both miscible and immiscible 
behavior. Immiscibility, where repulsion between species 
favours their spatial separation, has been observed [5- 
7, 14]. In recent years, many static and dynamical prop- 
erties of two-species BECs have been analysed. These 
include ground state structures [16-23], modulational 
instabilities [24-31], dark-bright solitons [32-43], vor- 
tices [11, 44-57], and the role of finite temperature [58- 
62]. In the limit of zero temperature, the mean-field of 
a single or two-component condensate is described by 
the Gross-Pitaevskii equation, in either single or cou- 
pled form, respectively. For immiscible two-component 
condensates under cylindrically symmetric trapping, the 
mean-field ground state has been shown to exist in a 
phase separated structure [16, 17, 63] where one compo- 
nent lies at the trap centre with the other lying at the 
periphery. This symmetry can be broken to give rise to 
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two separated side-by-side condensates [20, 23, 64, 65]. 

A recent two- species experiment with 87 Rb- 133 Cs con- 
densate mixtures (undertaken by some of the present au- 
thors) [5] revealed distinct regimes of density distribu- 
tions, depending on the relative numbers of 87 Rb and 
133 Cs atoms. This experiment exploits efficient sympa- 
thetic cooling of 133 Cs via elastic collisions with evap- 
oratively cooled 87 Rb atoms, initially in a magnetic 
quadrupole trap and subsequently in a levitated crossed 
dipole trap [66] . The large interspecies background scat- 
tering length of ~ 650 ao enables efficient sympathetic 
cooling but also gives rise to large inelastic three body 
losses [67]. This presents an obstacle to condensation 
at high densities which is overcome by fast evaporative 
cooling; this is achieved by reducing the dipole trap beam 
powers followed by tilting the trap using an applied mag- 
netic field gradient. Dual species condensates are pro- 
duced in the same trapping potential containing up to 
~ 2 x 10 4 atoms of each species. The optical dipole trap 
exerts harmonic trapping on the condensates. The mag- 
netic tilt contributes an additional weak linear potential 
in the tilt direction, which differs slightly between the 
species. Further trap effects, such as differential grav- 
itational sag, are also present. The combined result is 
that the potentials experienced by the two species have 
a slight offset in space. A dramatic spatial separation 
reveals this mixture to be immiscible as repulsive in- 
terspecies interactions dominate intraspecies interactions 
at the magnetic bias field of 22.4 G used in the experi- 
ment. As shown in Fig. 1 (a)-(b) the two-species BEC 
always forms one of three structures correlated to the 
atom number present in each condensate. Typical axial 
density profiles from each region (after time of flight ex- 
pansion) are shown in Fig. 1 (b). For Regions I and III, 
one of two possible symmetric cases is obtained: the Rb 
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FIG. 1. (Color online) (a) Experimental data for a quantum 
degenerate 87 Rb- 133 Cs mixture (reproduced from data in [5]). 
Depending on atom numbers, three distinct structures are 
observed represented here through triangles, squares and cir- 
cles (Regions I, II and III), (b) Experimental integrated axial 
density profiles corresponding to the filled symbols in (a), ob- 
served after time-of-flight expansion [68] and rescaled to the 
optical depth (OD) maximum, (c) Numerically calculated 
cylindrically symmetric ground state density profiles corre- 
sponding to the atom numbers for each of the filled points in 
(a) and the experimental density profiles shown in (b). (Solid) 
red curve — Rb; (dashed) blue curve — Cs. 



sits in the centre for Region I while the Rb is spatially 
split by the Cs in Region III. In Region II, the conden- 
sates adopt asymmetric density profiles, sitting side-by- 
side along the weaker axial direction of the trap. As we 
will see, these experimental profiles do not completely 
match the cylindrically-symmetric equilibrium solutions 
(presented in Fig. 1(c)). Furthermore, the rapid sympa- 
thetic cooling of this experimental system may lead to a 
situation where growth plays a determining factor in the 
final density structures formed. The current theoretical 
investigation is motivated by these considerations. 

In this paper, we study the equilibrium density struc- 
tures that arise in a 8T Rb- 133 Cs (referred to as Rb and 
Cs hereinafter) two-species condensate under harmonic 
trapping, with and without an asymmetrical perturba- 
tion between the two species (imposed via the addition 
of linear potentials) . Using zero temperature mean-field 
theory we characterize the crucial role of axial and trans- 



verse shifts in the harmonic traps of the two species in 
establishing such phase separation at equilibrium. 

In Sec. II we present our mean- field coupled Gross- 
Pitaevskii equations description of the two-species sys- 
tem, describe further the experimental parameters and 
the nature of the trapping potential used. In Sec. Ill A 
we present our ground state solutions of the coupled 
Gross-Pitaevskii equations under harmonic potentials. 
Supplementary results presented in Appendix A and Ap- 
pendix B explore the sensitivity of the numerical results 
to initial conditions and perform a dimensional analysis 
of the coupled equations, respectively. Sections III B and 
III C respectively examine the effect of adding to the un- 
derlying harmonic potentials further linear potentials in 
either the axial or one transverse direction. A combina- 
tion of both is studied for a set atom number in Sec. Ill D 
to see how the ground state density profiles of the system 
are influenced, with the gradients tailored to match ex- 
perimental results in Sec. HIE. The analysis presented 
here is based on the assumption that the profiles observed 
in the experiment are representative of the true equilib- 
rium state, something that we briefly comment upon in 
our concluding discussion section at the end of this paper 
(Sec. IV). 



II. THEORETICAL MODEL 

A. Coupled mean-field equations 

The Gross-Pitaevskii equation (GPE) is a nonlinear 
Schrodinger equation that describes the mean-field of a 
dilute Bose-Einstein condensate, first introduced to study 
vortex lines in an imperfect Bose gas [69-72]. Numerical 
simulations of the GPE, the simplest zero-temperature 
approach, frequently give excellent agreement with ex- 
periments [73-77]. In the case of a two-species BEC 
at absolute zero, one can use a set of coupled Gross- 
Pitaevskii equations (CGPEs) which should accurately 
describe equilibrium profiles in the limit of near-zero 
temperatures. These can be written as [16] 
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where ipi (x, y, z) and ip2 (x, z) are the mean-field wave 
functions for each condensate. Each wave function is 
normalized to its number of atoms independently i.e. 
J \ipi\ 2 dx dy dz = Ni(i — 1,2). The atomic mass 
and chemical potential for the zth species are respec- 
tively denoted by mi and The intraspecies interaction 
strengths are given by gu = AitTi 2 an/mi where an and 
<222 are the s-wave scattering lengths of each species re- 
spectively and the intercomponent coupling strength is 



3 



gi2 = 27rft 2 ai2/Mi2 where a\i is the interspecies scatter- 
ing length and Mi 2 the reduced mass [78]. The conden- 
sates are trapped in harmonic potentials of the form [76] 

Vi 0, y, z) = ^rrii (u 2 x{i) x 2 + ^ (i) y 2 + u 2 z{i) z 2 ^j , (3) 

where a;^)? ^(i) are the transverse trapping frequencies 
and uj z (j) the axial trapping frequency. This dynami- 
cal system has a substantial total parameter space, with 
even a restriction to cylindrical symmetry leaving eight 
dimensionless parameters that can in principle be inde- 
pendently varied (see Appendix B). Many of these pa- 
rameters will, in practice, be fixed in any given experi- 
mental configuration. Hence, for example, in the experi- 
mental configuration described in [5] , it is not possible for 
the distributions of Rb and Cs to be simply exchanged 
by changing the particle numbers (the most easily acces- 
sible handle to change the system's location in parameter 
space). This means, for example that the disagreement of 
Figs. 1(b) (i) and 1(c) (i) is unlikely to be due to incorrect 
atom counting. 

As discussed in [17], the two components can either 
overlap (miscible) or phase-separate (immiscible) de- 
pending on the relative strength of interactions between 
the two species. For an immiscible mixture, the strength 
of interactions for a homogeneous system requires [16, 77] 

g\ 2 > 9H922- (4) 

In a recent paper [79] it was shown that in inhomogeneous 
systems, phase separation can be suppressed even if this 
condition is fulfilled when the kinetic energy becomes im- 
portant. However, the system we consider strongly satis- 
fies the above criteria and lies deep within the immiscible 
regime. 

We use numerical simulations to obtain the three di- 
mensional ground states for the BEC mixture by solv- 
ing the CGPEs using the method of steepest descent [80] 
which amounts to simultaneously propagating (1) and (2) 
in imaginary time. As the initial trial solution we as- 
sume the absence of interspecies coupling and use the 
well-known Thomas-Fermi density profile for each con- 
densate rii = (fii — Vi)/gu [77]. The Thomas-Fermi pro- 
file becomes a good approximation for the exact den- 
sity profile when the kinetic energy contribution in the 
ground state is small. This occurs when Niaa/li ^> 1 
— where U = ^Ti/rriiUJi is the harmonic oscillator length 

of species i and Oi = {(^x^^y^^zii)) 1 ^ — a condition 
well-satisfied for our system. 

An added complexity of this mean-field model of two- 
species condensates is the occurrence, for certain param- 
eter regimes, of metastable steady state solutions of the 
system, which can be very close in energy to the true 
ground state. These solutions arise from different con- 
figurations of the two density profiles. The steady state 
solution obtained by imaginary time propagation is de- 
pendent on the initial state employed. We find that the 
Thomas-Fermi initial states detailed above consistently 



lead to the lowest energy solution, i.e. the ground state. 
We demonstrate the existence and behaviour of these 
metastable solutions in Appendix A, and discuss their 
presence in relation to our overall results in Sec. IV. 
However, it is important to note that all other results in 
this work relate to the true ground state of the system. 

Throughout the paper we present our obtained density 
profiles as either one-dimensional (ID) density profiles 
in the axial (z) direction, where the density has been 
integrated across both transverse directions, or as two- 
dimensional (2D) density profiles in the x-z plane, where 
the density has been integrated in the y direction. 

B. System parameters and test cases 

Following the experiment of Ref. [5] we consider a Rb- 
Cs BEC, denoting Rb as species 1 and Cs as species 2. 
The trap frequencies are taken to be cj^Rb) = w y (Rb) = 
2tt x 32.2 Hz, gj^Cs) = ^(Cs) = 27r x 40.2 Hz in the 
transverse directions and co^(Rb) = 2ttx 3.89 Hz, uj z (Cs) — 
2ir x 4.55 Hz in the axial direction. The intraspecies 
and interspecies scattering lengths are taken to be a R b = 
100 a [81], a Cs = 280 a [82] and a RbC s = 650 a [83]. 

Throughout this work we adopt harmonic oscillator 
units where time, length and energy are expressed in 
units of l/cD R b = 10 ms, Z R b = ^ n / raiCJ R b ~ 0.54 /urn 
and ftu) R b, respectively. 

The experiment of Ref. [5] observed three regimes of 
density structure, depending on the atom number in each 
species. These three regimes correspond to the triangles 
(Region I), squares (Region II) and circles (Region III) 
in Fig. 1 (a). We focus on one representative set of atom 
numbers from each structural regime: 

(i) ^Rb = 840 and TVcs = 8570 

(ii) N Rh = 3680 and TVcs = 8510 

(hi) N Rh = 15100 and TVcs = 6470 

These particular test cases are indicated by the three 
filled symbols in Fig. 1 (a) and correspond to the exper- 
imental images presented in Fig. 1 (b). 

C. Shifted trapping potentials 

In the experiment of Ref. [5] a magnetic tilt is applied 
to the otherwise harmonic potential to enhance evapora- 
tive cooling. This tilt is applied in one of the transverse 
directions, and results in a shift in relative trap centres 
by up to 3 microns. Additionally, the small difference 
in magnetic moment-to-mass ratio for each species, cou- 
pled with minute unavoidable misalignments of the dipole 
trap beams with respect to the magnetic potential, may 
result in offsets between the trap centres of up to 2 /im 
in all directions. To incorporate these shifts we add lin- 
ear potentials in the axial and one transverse direction 
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to species 1 such that its trapping potential takes on the 
modified form 

Vi = ^rai (ul {1) x 2 + Uy {1) y 2 + c^ (1) z 2 ) + ftxi + a z z. 

(5) 

The gradients for the linear potentials are given by a x 
and a z in the transverse and axial directions respectively. 
The corresponding expressions for the distance between 
the trap minima are 

&z — °% — and S x = — °^ — (6) 
rniu 2 z{1) ^i^(i) 

in the axial and transverse direction respectively. Our 
initial analysis focuses on symmetric ground state density 
profiles, obtained in the limit a x = a z = 0, for which the 
traps of each species are co-centered. 

III. EQUILIBRIUM PROFILES 

A. Symmetric System 

We first analyse equilibrium profiles in the absence of 
linear potentials, i.e. the trapping potentials for both 
species are centred relative to each other and are per- 
fectly harmonic. In such a case without any asymmetry 
in the trapping potential, we expect the obtained solu- 
tions to remain fully symmetric in all directions. Figure 1 
(c) shows the corresponding numerically-obtained inte- 
grated axial density ground states profiles. Our obtained 
states are all symmetric in space and phase-separated, as 
expected. For all three cases considered we observe the 
same qualitative structure in that the Cs cloud resides at 
the trap centre, with the Rb cloud surrounding it. This 
is the so-called 'ball-and-sheh" formation. Note that the 
observed preference in our numerical results for Cs to be 
centrally positioned is consistent with previous theoreti- 
cal studies where the tendency is for the component with 
the higher atomic mass to reside centrally [17, 84-86]. 
As we move from Region I to Region III, the qualitative 
structure does not change; only the relative amplitude of 
the condensates changes as increases and Nq s de- 
creases. 

Our results agree qualitatively with the experimental 
observations only for Region III, but not those obtained 
in Regions I and II. The experimental images include the 
presence of a broad thermal cloud, which is often large, 
and have undergone time-of-flight expansion, and so our 
comparison of density profiles is limited to the qualitative 
structural form only. In Region I, the experimental pro- 
files have the reverse structure to our numerical solutions, 
i.e. a central Cs condensate surrounded by Rb, whereas 
the experimental images for Region II are asymmetric in 
z. 

For the reasons discussed in Section II C, small shifts 
are present between the trapping potentials of the two 
species which introduce asymmetry into the system. We 




FIG. 2. (Color online) Ground states under a linear axial 
potential a z — 0.02 (ftu>Rb/fab)- Integrated (a) axial density 
profiles and (b) 2D density profiles, for the three cases corre- 
sponding to the filled symbols in Fig. 1 (a). (Solid) red curve 
— Rb; (dashed) blue curve — Cs. 

will now begin to investigate whether the introduction 
of these shifts and asymmetries may dictate the density 
structures that form and may enable us to reconcile the 
differences with the experimental results in Regions I and 
II. 



B. Axial Linear Potential 

We start by looking at how the addition of a linear po- 
tential in the axial direction influences the ground state 
solutions. Results presented here are based on a linear 
potential of gradient a z = 0.02 (fajRb/feb)- This cor- 
responds to an axial shift of S z = 0.9 /am in the trap 
centres, which is well within the experimental bounds 
detailed in Section II C. The value of the linear potential 
gradient employed here, and the value of the transverse 
linear potential employed in the next section, are chosen 
such that their combined result, presented in Sec. HIE, 
gives the best agreement to experiments that mean field 
theory can yield, while remaining within the bounds of 
the experimental uncertainties. 

The numerically-obtained ground state solutions are 
presented in Fig. 2 as both integrated (a) axial density 
profiles and (b) 2D density profiles. First consider the 
profile for Region III (Fig. 2 (hi)). The same structure 
remains from the symmetric system, i.e. central Cs sur- 
rounded by Rb, albeit now skewed slightly due to the 
linear potential. The profile for Region II (Fig. 2 (hi)) 
now jumps to an asymmetric side-by-side state, in qual- 
itative agreement with the experimental profile for this 
region (Fig. l(b)(ii)). 

For Region I (Fig. 2(i)) the density profile also be- 
comes asymmetric under the axial linear potential, and 
as such it remains inconsistent with the corresponding 
experimental observation (Fig. l(b)(i)). 

Increasing the linear potential past a critical value of 
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a z ~ 0.1 (^Rb/iRb) gives axially asymmetric density 
profiles for all three regimes while decreasing the gradi- 
ent below another critical value of a z c± 0.01 (h& R b/lRb) 
leads to the ball-in-shell structure for all three cases, i.e. 
central Cs surrounded by Rb. 



C. Transverse Linear Potential 

Next we look at the influence of an additional linear 
potential in one transverse direction, taken to be the x 
direction (no axial linear potential is applied here). We 
focus on a gradient of a x = 1.5 (hu R b/l R b) which corre- 
sponds to an axial shift of S x = 1.0 fim in the trap centres. 
This is comfortably within the experimental bounds de- 
tailed in Section II C. Figure 3 shows the corresponding 
integrated ground state density profiles. 

The integrated axial density profiles for Region III 
(Fig. 3 (a) (hi)) is again a qualitative match to that ob- 
tained experimentally. However, the actual 3D structure 
is now such that the Rb cloud lies offset transversely to 
the central Cs cloud, curving around it in the positive 
x half-plane, as visible in the corresponding integrated 
2D density plot (Fig. 3(b) (hi)). This is a subtly different 
structure to that observed in the previous sections where 
the Rb cloud was split either side of the Cs cloud. 

The ID density profile for Region II (Fig. 3(a)(ii)) 
shows similar results to the symmetric case where the Cs 
sits in the centre with Rb split axially into two distinct 
regions (and thus not consistent with the correspond- 
ing experimental profile). As visible in the 2D density 
plot (Fig. 3(b)(ii)), the outer Rb clouds become skewed 
towards the positive x direction, but not sufficiently to 
become linked to one side of the Cs cloud forming one 
Rb cloud. 

Importantly, for Region I, the integrated axial den- 
sity profile (Fig. 3(a) (i)) has undergone a marked change, 
with both condensates apparently overlapping. This does 
not contradict the presence of phase separation; as seen 
in the corresponding 2D density profiles (Fig. 3(b) (i)), 
the condensates phase separate transversely due to the 
transverse linear potential. Although this result is still 
different from the corresponding experimental profile 
(Fig. l(b)(i)), it is somewhat closer in that the Rb be- 
comes positioned in the centre in the axial direction. 

Increasing the linear potential past a critical value of 
a x ~ 3 (TicjRb/'Rb) gives transverse asymmetric density 
profiles for all three regimes while decreasing the gradient 
to values below approximately a x ~ 0.4 (hu R b/l R b) gives 
rise to a split in the Rb to surround the Cs positioned in 
the centre. 

So far we have presented numerical results for sym- 
metric trapping potentials and linear potentials in either 
the axial direction or a transverse direction. Results ob- 
tained for Region III qualitatively match those obtained 
experimentally. For Region II, we obtain a reasonable 
match through use of an additional axial linear potential. 
Through use of an additional transverse linear potential 
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FIG. 3. (Color online) Ground states under a linear transverse 
potential a x = 1.5 (ftu>Rb/fab)- Integrated (a) axial density 
profiles and (b) 2D density profiles, for the three cases corre- 
sponding to the filled symbols in Fig. 1 (a). (Solid) red curve 
— Rb; (dashed) blue curve — Cs. 



we have captured some, but not all, of the features from 
Region I. Given the difficulties in reproducing Region I, 
we will next focus on this case only, and explore the ex- 
tent to which a combination of both axial and transverse 
shifts, as relevant experimentally, enables a closer match 
to the experimental profiles. 



D. Combined Result: Axial and Transverse Linear 
Potential 

A summary of how different gradients of linear po- 
tentials in the axial and transverse directions affect the 
ground state density profiles is shown in Fig. 4. We dis- 
tinguish three distinct structures: a three peak profile 
(in which the Cs remains in the centre while the Rb is 
split axially into two distinct regions) , an axially side- by- 
side structure, or a transversely side-by-side structure. 
Starting from the symmetric three peaked distribution 
(bottom left in Fig. 4), we see that a small increase in 
a z gives rise to axially side-by-side density profiles: for 
a weak axial linear potential a z = 0.01 (hu R b/l R b), the 
ground state has switched to the axially side-by-side for- 
mation. On the other hand, a x needs to be increased 
more drastically to observe the switch of the ground 
state to a transversely side-by-side formation. This is 
due to competing effect of the condensate repulsion and 
the strong transverse trapping. When combining linear 
potentials in both axial and transverse directions, we see 
that for small a z < 0.04 (hu R b/l R b) the condensates 
favour a transverse side-by-side formation while this be- 
comes an axial side-by-side formation for larger values. 
For a z = 0.04 (huj Rh /l Rh ) and a x = 2.25 (hO Rh / l Rh ) , we 
see a combination of these side-by-side structures in that 
the Rb cloud lies diagonally to the side of the Cs cloud. 

For intermediate values of a x and a z , e.g. a x = 
0.4 (hO R b/lRb) and a z — 0.01 (ftc^Rb/^Rb), we see a trans- 
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FIG. 4. (Color online) Integrated 2D density profiles of the ground state as a function of the linear potential gradient in the 
axial (horizontal) and one transverse (vertical) direction applied to species l(Rb). Atom numbers correspond to point (i) in 
Fig. 1 (a). Red (black) — Rb; Blue (grey) — Cs. 



versely side-by-side formation in which the Rb is partic- 
ularly narrow and highly peaked, and the Cs features a 
small central density dip. We will see next that these 
profiles now give good agreement with the experimental 
observations for Region I. 



E. Matching of Experimental Profiles: Best 
Results of Mean Field Theory 

We now combine our previous analyses to demonstrate 
the extent to which addition of linear potentials in both 
the axial and one transverse direction matches the ex- 
perimental profiles for all three regimes for a given set of 
linear potentials. From Fig. 4, we can see that we require 
a z = 0.01 (huj Rh /l Rh ) or 0.02 (hu Rh /l Rh ) and a x = 0.4, 
1.5 or 2.25 (fttjRb/'Rb) to match the results for Region 
I. After a similar analysis for Region II and Region III, 
the experimentally relevant gradients of the linear po- 
tentials required to match all three simultaneously are 
a z = 0.02 (huj Rh /l Rh ) and a x = 1.5 (huj Rh /l Rh ) (hence 
justifying the use of those values in Sees. IIIB and III C). 
These linear potentials correspond to a displacement in 
the trap minima of S z = 0.9 /im and S x = 1.0 /im, both 
of which are well within experimental bounds detailed 
in Section II C. The numerical results for these gradi- 
ents are shown in Fig. 5 as both integrated axial den- 
sity profiles (Fig. 5(a)) and 3D isosurface plots of density 
(Fig. 5(b)). Our ground state results are now in very 
good qualitative agreement with the experimental pro- 



files (Fig. 1(b)) in all three regimes. For Region III, 
the ground state has the Rb cloud divided into two parts 
positioned either side of the central Cs cloud. For Re- 
gion II, we obtain the side-by-side formation in which 
the Cs cloud sits to the right of Rb cloud. For Region 
I, the ground state features the Rb cloud to be centrally 
located in the axial direction, but shifted transversely, 
while the Cs has a small density dip at the centre. 

The central density dip in the Cs profile for Region 
I is more pronounced in the experimental observations, 
e.g. Fig. 1(b) (i), than in our above results. An inherent 
feature of solving the coupled Gross-Pitaevskii equations 
for an immiscible two-species BEC is a sensitivity to the 
initial trial wavefunction. All of our results presented so 
far have been based on TF initial trial wavefunctions (for 
the quoted condensate atom numbers for each species), 
as described in Section II A. By their nature, the TF 
profiles tend to be broadly distributed in space, and this 
tends to favour a more broad distribution of condensates 
in the final static solution obtained. We find that em- 
ploying an initial distribution for the Rb cloud which 
is tightly localized at the origin yields static solutions 
which feature a localized Rb cloud and a slightly more 
prominent density dip in the adjacent Cs cloud, in closer 
agreement with the experimental profiles for Region I. 
These numerical results are presented in Appendix A. 

We have also looked at introducing axial asymmetry 
(without permanent trap shifts) through shifts in our 
numerical initial conditions. Specifically, the TF initial 
conditions for each species are initially offset along the 
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FIG. 5. (Color online) Ground states under an axial linear 
potential a z = 0.02 (/icjRb/^Rb) and a transverse linear po- 
tential a x = 1.5 (ftcjRb/feb). (a) ID density profiles and (b) 
3D isosurface plots, each corresponding to the filled symbols 
in Fig. 1 (a). (Solid) red curve — Rb; (dashed) blue curve — 
Cs. 

z-axis. Similarly to the use of the linear potential, this 
initial offset could be tailored to reproduce the experi- 
mental results to a similar degree of accuracy. 

We have additionally simulated the expansion of the 
static solutions following the sudden removal of all trap- 
ping potentials. This expansion is performed experimen- 
tally prior to imaging. Our analysis showed that expan- 
sion does not affect the structures formed. The overall 
phase separation features appeared to be captured very 
well under the assumption made here that the profiles ob- 
served in the experiments are the true equilibrium profiles 
and that these profiles are dominated by their respec- 
tive condensate component, with thermal clouds simply 
modifying these profiles by the addition of characteristic 
thermal tails. 



IV. DISCUSSION 

We investigated the 8T Rb- 133 Cs ground state density 
profiles corresponding to the parameters of a recent ex- 
periment [D. J. McCarron et al, Phys. Rev. A 84, 
011603 (2011)] by means of the simplest possible zero- 
temperature mean-field theory consisting of two coupled 
Gross-Pit aevskii equations. Density profiles obtained in 
perfectly symmetric traps were found not to match the 
experimental results. Analysing the experiment more 
carefully, we proceeded to add weak linear potentials in 
the axial and one transverse direction accounting for an- 
ticipated experimental offsets (of around 1 fim) in the 
trap centres for the two species. Even weak linear poten- 
tials can give rise to dramatically different density pro- 
files. Importantly, this allows us to obtain the observed 
asymmetric profiles. By tailoring the gradients of the 
linear potentials, we find our simulations to qualitatively 
match structural regimes seen experimentally, when fo- 



cusing on condensate separation features and overlooking 
the appearance of thermal tails. 

The analysis presented in this work was based on equi- 
librium density profiles. While we demonstrated good 
overall agreement with the experimentally-reported pro- 
files, we also found that a change in the initial condi- 
tions of the simulations, e.g. one of the components be- 
ing more tightly localised in the centre, could affect the 
final equilibrated profiles, as numerous metastable states 
(of comparable, but not identical, energies) exist for each 
configuration. Such a situation could for example arise 
in the early stages of coupled growth under some param- 
eter regimes. In the experiments, as the two species were 
sympathetically cooled, the initial number of condensate 
atoms within each species (or the sequence by which 
growth proceeded) was not accurately known. Moreover, 
the density profiles were typically measured after a vari- 
able hold time, without necessarily guaranteeing that the 
structures observed were indeed true equilibrium states 
(as opposed to some long-lived metastable steady- states) , 
for which a detailed analysis of growth dynamics would 
have been required. Preliminary investigation of cou- 
pled Gross-Pitaevskii equations with phenomenological 
damping undertaken by us indeed revealed different fea- 
tures during growth, depending on both initial conditions 
and growth parameters; more importantly, however, such 
simulations showed that after sufficient evolution time, 
the condensate in one or the other species disappeared, 
a feature which is in qualitative agreement with the ex- 
periments, which detected only one condensate (of ei- 
ther species) in some measurements. The study of cou- 
pled two-component condensate growth is an interesting 
topic that will be studied in more detail in future work. 
Similar non-equilibrium conclusions have been reached 
by another group [87] using such equations additionally 
modified by the presence of stochastic noise mimicking 
thermal fluctuations, which additionally allows for the 
appearance of spontaneous structures during growth. 
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Appendix A: Sensitivity to Initial Conditions 

Instabilities often arise when solving equations numer- 
ically. Here, we show how these can arise during conver- 
gence of the CGPE (during imaginary time propagation) 
for various initial conditions. For simplicity we focus on a 
ID two-species system. The relation between the ID and 
3D mean- field equations are given by Uu = <fo/27rZ?j_, 
U12 = #12 A + Vi&D) = fa - hWi( Xi y) and 

li± = y/h/miWi( Xj y) is the transverse harmonic oscilla- 
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FIG. 6. (color online) Integrated axial ground state den- 
sity profiles, with N = Ni = N2, uji =002, mi — rri2 and 
U22 = l.OlC/n, U12 = 1.52C/n, such that the immiscibility 
criterion is satisfied. Columns correspond to Gaussian, TF 
and homogeneous initial conditions for imaginary time prop- 
agation, (a) Un = 6, N = 150. (b) Uu = 6, N = 200. (c) 
Uu = 6, N = 2000. (d) Uu = 1, N = 200. Solid blue curve 
- species 1; Dashed red curve - species 2. 
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FIG. 7. (Color online) Same parameters for uji, rrii, U22 and 
U 12 as in Fig. 6 and N = 2000, Uu = 6. (a) Number of 
interfaces between species plotted against the width of the 
gaussian trial solution £. (b) Numerical initial conditions for 
points (i)(C = 1) and (ii)(C = 10) respectively, (c) Ground 
state density profiles for points (i) and (ii) respectively. Solid 
blue curve - species 1; Dashed red curve - species 2. 



tor length. Similarly to [20], we take N = Ni = N 2 , 
uji = UJ2 1 mi = 1712 and U22 = l.OlC/n, U12 = 1.52/7n. 

We consider the two-species condensates to be con- 
fined by co-centered harmonic traps. We thus antici- 
pate symmetric density profiles to emerge. We perform 
imaginary time propagation of the CGPEs subject to dif- 



ferent initial conditions. Convergence is decided when 
the fractional difference in energy between time steps is 
lower than 10 -9 . We consider (i) the Gaussian ground 
state harmonic oscillator solution for each species, (ii) a 
Thomas- Fermi (TF) initial state for each species, and (hi) 
a homogeneous initial state (uniform density) for each 
state. Each initial condition is suitably normalized to 
the desired N. 

In Fig. 6, we compare the steady state density pro- 
files that emerge through each of these initial conditions. 
The first row, Fig. 6(a), corresponds to the parameters 
U u — 6 an d N = 150. The oscillator length I is given by 
The Gaussian initial condition gives rise to 
a steady state which features many corrugations in den- 
sity. The TF initial state gives rise to a central cloud of 
species 1 surrounded by species 2. This is consistent with 
previously obtained solutions [20]. For the homogeneous 
initial state, we obtain a steady state which is similar 
to the TF result but features additional small clouds of 
species 1 at the periphery. 

The fractional difference in energies for each of these 
systems shows that the ground states given by the TF 
initial conditions have the smallest energies followed by 
the homogeneous initial condition and finally the Gaus- 
sian. We thus identify the state derived from the TF 
initial condition to be the ground state of the system. 
The difference is energies, however, is small and typi- 
cally less than 10%. The additional corrugations arising 
from the Gaussian and homogeneous initial conditions 
are attributed to a modulational instability of the con- 
densates early in the imaginary-time propagation. This 
leads to the formation of density corrugations in the sys- 
tem. Imaginary time propagation then converges to one 
of several metastable states with these corrugated forms. 

For a slightly increased atom number, A/", Fig. 6 (b), 
we find that the TF-based solution becomes identical to 
the homogeneous-based solution, both featuring small 
clouds of species 1 at the periphery. The Gaussian- 
based solution is essentially unchanged. For a signifi- 
cantly larger atom number, Fig. 6 (c), all profiles become 
wider due to the increased repulsion in the system. The 
TF and homogeneous-based solutions maintain the same 
structural form as above, albeit broader. However, the 
Gaussian-based solution features a much greater number 
of density corrugations. Finally, in Fig. 6 (d), we plot the 
results for when the interaction strengths between atoms 
are reduced compared to Fig. 6 (b). In this regime of 
weak interactions, all three initial conditions lead to the 
same solution in which a cloud of species 1 is surrounded 
by species 2. 

The Gaussian-based solution features the greatest sen- 
sitivity. In order to parameterize these effects further we 
introduce the number of interfaces as a measure of the 
number of corrugations in the system. We define our 
Gaussian initial condition (in dimensionless form) as, 

|^|2 =e -^/2C 2 . (A1) 

Above, we employ the ground harmonic oscillator solu- 
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tion corresponding to £ = 1. Here we will vary the Gaus- 
sian width ( to assess its role on the steady state solution 
obtained. In Fig. 7(a) we plot the number of interfaces 
in the steady state solution as a function of £. As £ 
is increased (the Gaussian initial state becomes wider), 
the number of interfaces decreases exponentially until the 
minimum is reached, matching the number of interfaces 
in the ground state (which corresponds to the TF and 
homogeneous-based solutions in Fig. 6(c)). In Fig. 7(b) 
and (c), we show the initial conditions and ground state 
density profiles for two cases: (i) £ = 1 and (ii) ( = 10 
respectively. These results show firstly that a range of 
metastable states are possible in the system, character- 
ized by an increased number of density corrugations over 
the ground state. Secondly, the results suggest that the 
size of the initial trial solution relative the ground state is 
key to determining which state is obtained via imaginary 
time propagation. 

In 3D the solutions obtained are also sensitive to the 
initial conditions, although considerably less so than in 
ID. Consider the test case in Region I and fully symmet- 
ric trapping (a x = a z =0). Using TF profiles as initial 
states we obtained the integrated axial density profile 
shown in Fig. 1(c) (i) in which the Rb sits either side of 
the central Cs cloud. By beginning instead with a very 
narrow Gaussian profile for the Rb while retaining the TF 
profile for Cs (assuming essentially here that the Cs con- 
densate grew first) we can obtain a metastable solution 
whose integrated axial profile features the Rb sitting in at 
the trap centre and a small density dip in the ambient Cs 
cloud, in qualitative agreement with the corresponding 
experimental profile. These solutions have higher energy 
than that of the solutions observed in Fig. 1 (c) (i) and 
so therefore are not the 'true' ground states. However, 
in the presence of trap shifts (a x , a z ^ 0), we essentially 
regain the true solution presented previously in Figure 
5(a) (i), i.e. the effect of the initial numerical state be- 
comes washed out. Given we that achieve essentially no 
improvement in our final comparison to the experimental 
profiles by changing the initial states, we do not present 
any further results on this. 



Appendix B: Dimensionless Analysis of Coupled 
Gross— Pitaevskii Equations. 



2m 



-V 2 + ^[A 2 (* 2 + ^) + , 2 ] 

+011 l^l| 2 +012 |^2| 2 } ^1 



dt 



+022 |^2| 2 +012 |^l| 2 }^2, 



where ga and gi2 are given in Sec. II A, and Ai and A2 
are the trap aspect ratios of component 1 and 2, respec- 
tively. Note also that the chemical potentials fii and [12 
are not relevant to the real-time dynamics of the CG- 
PEs (contributing only a global phase to the condensate 
wavefunctions), and are thus not explicitly stated. We 
choose a set of symmetrized harmonic units, codified as 
ft = = ^ojiuj2 = 1. This means we have time, 

length, and energy units of 

r =(cjicj 2 )" 1/2 , 



ft 



yraira 2 co>iCo> 2 



e =ft v / cjicj 2 - 

Rescaled renderings of the coupled GPEs (where we also 
normalize the wavefunctions ipi and ip2 to 1, in order to 
make the dependence on the particle numbers Ni and N2 
of the two species more explicit), expressed in terms of a 
minimal number of parameters, are then 



.0^1 f 7V 2 1 r 2 2 2 21 



+ 7«n |^i I 



7 2 + l 
2 7 




+ 



7 2 + 1\ OL\2 I / .2 , a 22 I 



2 7 



-|^l| 2 + — |^2| 2 ^2, 

V 7 



where 



In this work we focus on the experimental Rb-Cs sys- 
tem and as such we proceed with an analysis based on 
parameters for this system. However, more generally, the 
CGPEs can be reduced to a dimensionless form. Here we 
outline this reduction and highlight some key features 
that emerge. We start from the CGPEs, where the po- 
tential is provided by cylindrically symmetric harmonic 
trapping potentials with common minima, 



7 = 




AV 
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Using the total particle number N — N\ + N2 , the inter- 
action coefficients can alternatively be phrased as 

N 

an 



■ 2 + l' 
T]N 

J] 2 



which is a particularly natural description if the two 
species are simply different internal states of the same 
atom (in which case 7 = 1, simplifying the system of 
equations further). We therefore have eight independent 
dimensionless parameters (including Ai and A2). 
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